% 円制限三体問題に変更
function [f_xnext,f_dxn_dx,f_dxn_du,n,m,N,alphac,epsilon,baru,T,deltat,mu_e,mu_m,mu,r_lagrange_1, r_lagrange, r_lagrange_3] = make_problem()

alphac = 1;%任意設計パラメータ
epsilon=0.000025;%default epsilon = 5*10^(-3)
baru=8000;

n = 4;% 状態 x = [x, y, v_x, v_y]T の次元
m = 2;% 入力 u = [Ux, Uy]T の次元
x = sym('x',[n,1],'real');% シンボリック変数 'x' を作成する．シンボリック演算に用いるための値や数式を保存するための変数.
u = sym('u',[m,1],'real');% シンボリック変数 'u' を作成する．シンボリック演算では数式に対して数値ではなく数式のままの状態で行われる．
T  = 0.75;%0.75
N = 2000;%離散化数
deltat = T/N;

mu_e = 3.98600441*10^14;
mu_m = 4.90279898*10^12;
mu = mu_m/(mu_e+mu_m);
r_lagrange_1 = 0.150934272283726; %L1 lagrange point is 1-mu-r_lagrange_1
r_lagrange = 0.167832730863978; %L2 lagrange point is 1-mu+r_lagrange
r_lagrange_3 = 0.0070879373688714; %L3 lagrange point is -mu-(1-r_lagrange_3)

%% yamamoto runge kutta
F_c =[x(3);
      x(4);
      u(1) + 2*x(4) + x(1) - (1-mu)*(x(1)+mu)/((x(1)+mu)^2+x(2)^2)^(3/2) - mu*(x(1)-1+mu)/((x(1)-(1-mu))^2+x(2)^2)^(3/2) ;
      u(2) - 2*x(3) + x(2) - (1-mu)*x(2)/((x(1)+mu)^2+x(2)^2)^(3/2) - mu*x(2)/((x(1)-(1-mu))^2+x(2)^2)^(3/2) ];
F = matlabFunction(F_c,'Vars',{x,u});
% Runge-Kutta
k1 = @(x,u) F(x,u);
k2 = @(x,u) F(x + deltat*k1(x,u)/2,u);
k3 = @(x,u) F(x + deltat*k2(x,u)/2,u);
k4 = @(x,u) F(x + deltat*k3(x,u),u );
pre_f =( x + deltat*(k1(x,u) + 2*k2(x,u)+ 2*k3(x,u) + k4(x,u))/6);
f_xnext =@(x,u) x + deltat*(k1(x,u) + 2*k2(x,u)+ 2*k3(x,u) + k4(x,u))/6;% x(k+1)を計算してくれる．

pre_delf_delu = jacobian(pre_f,u);
pre_delf_delx = jacobian(pre_f,x);

f_dxn_dx = matlabFunction(pre_delf_delx,'Vars',{x,u});
f_dxn_du = matlabFunction(pre_delf_delu,'Vars',{x,u});

end